O.4 Frontiers OXT

OXT and Psichological variables - By Gender

1 Read and Merge

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Read 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
wd <- "C:/Users/Mihai/Desktop/O.4 prealabil pt Frontiers/O.4 Date Psiho&OXT"
setwd(wd)

Data_OXT <- rio::import(file.path(wd, "\\OXT", "O.4 OXT Data.xlsx"))

Data_psih <- rio::import(file.path(wd, "\\Psiho", "DATE O4Cl 16.11.2019 procesat.xlsx"), which = "26martie2019")
Data_psih <- Data_psih[, c(1:6, 8:13)]

Data_Neo <- rio::import(file.path(wd, "\\Psiho", "DATE O4Cl 16.11.2019 procesat.xlsx"), which = "Neo scored")
Data_BDI <- rio::import(file.path(wd, "\\Psiho", "DATE O4Cl 16.11.2019 procesat.xlsx"), which = "BDI scored")
Data_STAI <- rio::import(file.path(wd, "\\Psiho", "DATE O4Cl 16.11.2019 procesat.xlsx"), which = "STAI scored")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Clean & Merge
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Data_STAI <- 
  Data_STAI %>% 
  dplyr::filter(!is.na(ID)) %>%                                                                # no ID, no data
  dplyr::filter(rowSums(is.na(.)) < 2) %>%                                                     # filter out rows (no more than 1 NA on row)
  mutate(ID = stringr::str_replace(ID, "\\s", "|")) %>%                                        # replace first whitespace with | and use it to separate
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\|") %>%                              # separate on only first whitespace that was replaced
  dplyr::select(-Experim)

Data_BDI <- 
  Data_BDI %>% 
  dplyr::filter(!is.na(ID)) %>%                                                                # no ID, no data
  dplyr::filter(rowSums(is.na(.)) < 1) %>%                                                     # filter out rows (not 1 NA)
  mutate(ID = stringr::str_replace(ID, "\\s", "|")) %>%                                        # replace first whitespace with | and use it to separate
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\|") %>%                              # separate on only first whitespace that was replaced
  dplyr::select(-Experim)

Data_Neo <- 
  Data_Neo %>% 
  dplyr::filter(!is.na(ID)) %>%                                                                # no ID, no data
  dplyr::filter(rowSums(is.na(.)) < 6) %>%                                                     # filter out rows (no more than 5 NA on row)
  mutate(ID = stringr::str_replace(ID, "\\s", "|")) %>%                                        # replace first whitespace with | and use it to separate
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\|") %>%                              # separate on only first whitespace that was replaced
  dplyr::select(-Experim)

Data_psih <- 
  Data_psih %>% 
  dplyr::filter(!is.na(ID)) %>%                                                                # no ID, no data
  dplyr::filter(rowSums(is.na(.)) < 8) %>%                                                     # filter out rows (no more than 7 NA on row)
  tidyr::separate(Conditia,  c("Nr_zi", "Conditia"), "\\s+") %>%                               # split on white space
  mutate(ID = stringr::str_replace(ID, "\\s", "|")) %>%                                        # replace first whitespace with | and use it to separate
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\|") %>%                              # separate on only first whitespace that was replaced
  dplyr::rename("Nume" = Nume_Prenume) %>% 
  dplyr::select(-Experim)

Data_OXT <-   
  Data_OXT %>% 
  tidyr::separate(ID,  c("ID", "Ziua", "Nr_zi", "Proba"), "\\s+") %>%                          # split on white space
  select(ID, Nr_zi, Proba, OXT) %>%
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\/") %>%                              # separate on /
  dplyr::select(-Experim) %>% 
  dplyr::rename("PrePost" = Proba) %>% 
  mutate(PrePost = forcats::fct_recode(PrePost, "Pre" = "A", "Post" = "B")) %>%                # Proba A = Pre, B = Post on same day
  spread(key = PrePost, value = OXT) %>%
  dplyr::rename_at(vars(Pre, Post), ~ c("OXT_Pre", "OXT_Post"))


Data_merged <- dplyr::left_join(Data_psih, Data_OXT, by = c("ID", "Nr_zi"))                 # ID 40, 59 are incomplete in OXT -- dplyr::full_join


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Make Wide Dataframe of merged
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Check if ids have > 1 row of data (empty ID have only 1 row; thei have only 1 Condition) 
# Careful! This function modfies the datasets in the global envinronment
delete_empty_id <- function(df){
  list_empty_id <- 
    df %>%
    dplyr::group_by(ID) %>%
    dplyr::summarise(row_count = n()) %>%
    dplyr::rename("empty_id" = ID) %>%
    mutate(delete_id = if_else(row_count < 2, TRUE, FALSE)) %>%
    filter(delete_id == TRUE)
  
  df_modif <- 
    df %>%
    filter(!ID %in% list_empty_id$empty_id)
  
  if(!identical(df, df_modif)){
    df <- deparse(substitute(df))
    cat("Deleting from ", print(as.name(df))); print(list_empty_id)                    # print out which ids are deleted from which dataset
    assign(df, df_modif, envir = globalenv())                                          # assign modified df to original dataset from Global
  }else cat("No empty datasets. Nothing to delete")
}
delete_empty_id(Data_merged) 

Data_merged Deleting from Data_merged

# Dataframe for regular analyses
Data_merged_wide <-
  Data_merged %>%
  select(-c(Varsta, Gen)) %>%            # exclude Gen, Varsta: merge after because it interferes with spread()
  gather("variable", "value", c(Data, Nr_zi, Vas_Stres_Pre:OXT_Post), -c(1:2)) %>%       # Conditia needs to be outside
  unite("united_var", c(variable, Conditia), sep = "_") %>%
  spread(united_var, value) %>%
  mutate_at(vars(5:22), list(~as.numeric(as.character(.))))

# sum(is.na(Data_merged[, 6:15])) == sum(is.na(Data_merged_wide[, 5:22]))    # check if nr of NA is the same after melt


# Full Wide Dataframe  (HERE We APPLY EXCLUSION CRITERIA) 
Data_merged_wide_all <- 
  Data_merged[, c("ID", "Nume", "Varsta", "Gen")] %>%                        # need to add these back, but just the main unduplicated row
  dplyr::distinct(ID, .keep_all = TRUE) %>%
    dplyr::left_join(., Data_merged_wide, by = "ID") %>%
    dplyr::rename("Nume" = Nume.x) %>%
    dplyr::select(-Nume.y) %>% 
    dplyr::left_join(., Data_BDI, by = "ID") %>%
    dplyr::left_join(., Data_STAI, by = "ID") %>%
    dplyr::left_join(., Data_Neo, by = "ID") %>%
    mutate_at(vars(23:60), list(~as.numeric(as.character(.)))) %>%
    mutate(Varsta = as.numeric(as.character(Varsta))) %>% 
    filter(!ID %in% c(19, 43)) %>%           # 55 subs remain after this
    filter(Varsta < 40) %>%                  # 53 subs remain after this
    filter(ScorBDI < 30) %>%                 # 50 subs remain after this
    dplyr::arrange(ID)


# Full Long Dataframe
vars_OglEcran <- grep("OGL|ECRAN", colnames(Data_merged_wide_all))

Data_merged_long_all <-
  Data_merged_wide_all %>%
  gather(variable, value, vars_OglEcran, -c(ID)) %>%                                      # gather Data and Nr_zi as these are OglEcran level
  tidyr::separate(variable,  c("variable", "OglEcran"), "_(?=[^_]+$)") %>%                # split only on last "_"
  spread(variable, value) %>% 
  rename_at(vars(vars_OglEcran), list(~stringr::str_remove_all(., c("_OGL|_ECRAN")))) %>%          
  arrange(ID)

vars_PrePost <- grep("Pre|Post", colnames(Data_merged_long_all))

Data_merged_long_all <-
  Data_merged_long_all %>%
  pivot_longer(vars_PrePost, names_to = c(".value","PrePost"), names_pattern = "(.*)_(.*)") %>%     # magic     
  mutate_at(vars(48:50), list(~as.numeric(as.character(.)))) %>%
  arrange(ID)


# Females and Males Dataframes
Data_fem <- 
  Data_merged_wide_all %>%
  filter(Gen == "f")
Data_masc <- 
  Data_merged_wide_all %>%
  filter(Gen == "m")

# Difference Scores Dataframe
Data_difscores <-
  Data_merged_wide_all %>%
  mutate(Diff_OXT_OGL = OXT_Post_OGL - OXT_Pre_OGL,
         Diff_OXT_ECRAN = OXT_Post_ECRAN - OXT_Pre_ECRAN,
         Diff_Vas_Stres_OGL = Vas_Stres_Post_OGL - Vas_Stres_Pre_OGL,
         Diff_Vas_Stres_ECRAN = Vas_Stres_Post_ECRAN - Vas_Stres_Pre_ECRAN, 
         Diff_Vas_Bine_OGL = Vas_Bine_Post_OGL - Vas_Bine_Pre_OGL,
         Diff_Vas_Bine_ECRAN = Vas_Bine_Post_ECRAN - Vas_Bine_Pre_ECRAN,
         Diff_IOS_OGL = IOS_Post_OGL - IOS_Pre_OGL,
         Diff_IOS_ECRAN = IOS_Post_ECRAN - IOS_Pre_ECRAN
  ) 

1.1 Dataset - Long Format

Data_merged_long_all %>%
  dplyr::select(-Nume) %>%
    DT::datatable(
      extensions = 'Buttons',
      options = list(pageLength = 10,
                     scrollX='500px',
                     dom = 'Bfrtip',
                     buttons = c('excel', "csv")))

1.2 Dataset - Wide Format

Data_merged_wide_all %>%
  dplyr::select(-Nume) %>%
    DT::datatable(
      extensions = 'Buttons',
      options = list(pageLength = 10,
                     scrollX='500px',
                     dom = 'Bfrtip',
                     buttons = c('excel', "csv")))

2 Descriptives

2.1 Descriptives - Gen Varsta

describe(Data_merged_wide_all$Varsta)

prop.table(table(Data_merged_wide_all$Gen))

   f    m 
0.54 0.46 

t.test(Data_merged_wide_all$Varsta ~ Data_merged_wide_all$Gen) %>%
  tidy()

wilcox.test(Data_merged_wide_all$Varsta ~ Data_merged_wide_all$Gen) %>%
  tidy()

cannot compute exact p-value with ties

psych::describeBy(Data_merged_wide_all$Varsta, group = Data_merged_wide_all$Gen)

 Descriptive statistics by group 
group: f

------------------------------------------------------------------------------------------------------------ 
group: m

2.2 Plots

## Pie chart
Data_merged_wide_all  %>%
  mutate(Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "f", "masculin" = "m")) %>%
  group_by(Gen) %>%
  dplyr::summarise(counts = n()) %>%
  mutate(prop = round(counts*100/sum(counts), 1),
         lab.ypos = cumsum(prop) - .5*prop,
         Percent = paste0(prop, " %")) %>% 
  ggpubr::ggpie(x = "prop", label = "Percent",
                fill = "Gen", color = "white", 
                lab.pos = "in", lab.font = list(color = "white"),
                palette = "grey")

## Dodged Bar plot of Age and Gender
Data_merged_wide_all %>%
  mutate(Varsta = as.numeric(as.character(Varsta))) %>%
  mutate(Varta_categ = cut(Varsta, 
                           breaks=c(-Inf, 25, 30, 35, 40), 
                           labels=c("20-24","25-29","30-34", "35-39"), 
                           right = FALSE)) %>%  
  mutate(Varsta = as.factor(Varsta),
         Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "f", "masculin" = "m")) %>%
  dplyr::count(Varta_categ, Gen, .drop = FALSE) %>%         # Group by, then count number in each group (dont drop 0 counts)
  mutate(pct = prop.table(n)) %>%                           # Calculate percent within each var
    ggplot(aes(x = Varta_categ, y = pct, fill = Gen, label = scales::percent(pct))) + 
      geom_col(position = position_dodge(preserve = "single"), stat = "identity") +    # Don't drop zero count
      geom_text(position = position_dodge(width = .9),      # move to center of bars
                vjust = -0.5,                               # nudge above top of bar
                size = 3) + 
      scale_y_continuous(labels = scales::percent) +
      ggtitle("") +
      xlab("Varsta") + ylab("Percentage %") + 
      guides(fill = guide_legend(title = "Gen", ncol = 1)) + 
      scale_fill_grey(start = 0.8, end = 0.2, na.value = "red", aesthetics = "fill") +
      theme(legend.position = "right", legend.direction = "vertical", 
            legend.justification = c(0, 1), panel.border = element_rect(fill = NA, colour = "black"))

2.3 t test Varsta by Gen

# Gen
Data_merged_wide_all  %>%
  mutate(Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "f", "masculin" = "m")) %>%
  group_by(Gen) %>%
  dplyr::summarise(counts = n())

# varsta
Data_merged_wide_all %>%
  mutate(Gen = as.factor(as.character(Gen))) %>%
  mutate(Varsta = as.numeric(as.character(Varsta))) %>%
  tidystats::describe_data(Varsta, na.rm = TRUE)

# varsta by Gen
Data_merged_wide_all %>%
  mutate(Varsta = as.numeric(as.character(Varsta))) %>%
  dplyr::group_by(Gen) %>%
  tidystats::describe_data(Varsta, na.rm = TRUE) %>%
  knitr::kable(caption = "Varsta by Gen", format = "pandoc", digits = 2)

Varsta by Gen

var	Gen	missing	n	M	SD	SE	min	max	range	median	mode	skew	kurtosis
Varsta	f	0	27	25.44	5.49	1.06	20	37	17	24	20	0.78	2.30
Varsta	m	0	23	26.00	4.67	0.97	20	39	19	25	23	0.81	3.62

# t test varsta by Gen
Data_merged_wide_all %>%
  mutate(Varsta = as.numeric(as.character(Varsta))) %>%
  t.test(Varsta ~ Gen, data = .) %>% 
  tidy()

3 Define Functions

subchunkify <- function(g, fig_height=7, fig_width=5) {        # THIS FUNC CRASHES PANDOC CONVERSION WITH ERROR 99 - DONT USE HERE
  g_deparsed <- paste0(deparse(
    function() {g}
  ), collapse = '')
  
  sub_chunk <- paste0("
  `","``{r sub_chunk_", floor(runif(1) * 10000), ", fig.height=", fig_height, ", fig.width=", fig_width, ", echo=FALSE}",
  "\n(", 
    g_deparsed
    , ")()",
  "\n`","``
  ")
  
  cat(knitr::knit(text = knitr::knit_expand(text = sub_chunk), quiet = TRUE))
}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Define Function for mining correlations
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Function for p-value significance -- both for func_ancova_multibox(), Get_Top_Relationships() and Correlations_With_One()
stars_signif <- function(pval) {
  stars = "ns"
  if(pval <= 0.001)
    stars = "***"
  if(pval > 0.001 & pval <= 0.01)
    stars = "**"
  if(pval > 0.01 & pval <= 0.05)
    stars = "*"
  if(pval > 0.05 & pval <= 0.1)
    stars = "."
  stars
}

## Function that returns correlations of all variables in descending order.
# Arg for threshold with default at .3 will keep only correlantions above .3 and below -.3. Also has threshhold for p-value. 
Get_Top_Relationships <- function(data_set, 
                                  correlation_abs_threshold=0.3,
                                  pvalue_threshold=0.05) {
  require(psych)
  require(dplyr)
  feature_names <- names(data_set)
  # strip var names to index for pair-wise identification
  names(data_set) <- seq(1:ncol(data_set))
  # calculate correlation and significance numbers
  cor_data_df <- psych::corr.test(data_set)
  # apply var names to correlation matrix over index
  rownames(cor_data_df$r) <- feature_names
  colnames(cor_data_df$r) <- feature_names
  # top cor and sig
  relationships_set <- cor_data_df$ci[,c('r','p')]
  # apply var names to data over index pairs
  relationships_set$feature_1 <- feature_names[as.numeric(sapply(strsplit(rownames(relationships_set), "-"), `[`, 1))]
  relationships_set$feature_2 <- feature_names[as.numeric(
    sapply(strsplit(rownames(relationships_set), "-"), `[`, 2))]
  relationships_set <- dplyr::select(relationships_set, feature_1, feature_2, r, p) %>% dplyr::rename(correlation = r, p.value = p)
  # return only the most insteresting relationships
  return(filter(relationships_set, abs(correlation) > correlation_abs_threshold &
                  p.value < pvalue_threshold) %>% 
        arrange(p.value) %>%
        mutate(p.signif = sapply(p.value, function(x) stars_signif(x))))
}

## Function that returns all correlation between numeric variables and one specific variable
Correlations_With_One <- function(data_set,
                            variable,
                            correlation_abs_threshold=0.3,
                            pvalue_threshold=0.05) {
  require(psych)
  require(dplyr)
  # use all numeric columns only
  numeric_cols <- unlist(lapply(data_set, is.numeric))
  data_set <- data_set[, numeric_cols]                               
  # calculate correlation and significance numbers
  cor_data_df <- psych::corr.test(data_set[, names(data_set) != variable], data_set[, variable], minlength = 20, adjust="none")
  # top cor and sig
  relationships_set <- as.data.frame(cbind(cor_data_df$r, cor_data_df$p))     # same as  cor_data_df$ci[,c('r','p')]
  relationships_set <- tibble::rownames_to_column(relationships_set, "Variable")   # relationships_set$Variable <- rownames(relationships_set)
  colnames(relationships_set) <- c("Variable", "correlation", "p.value")
  # return only the most insteresting relationships
  cat("#### Correlations with ", variable, "\n")
  return(filter(relationships_set, abs(correlation) > correlation_abs_threshold &
                  p.value < pvalue_threshold) %>% 
           arrange(p.value) %>%
           mutate(p.signif = sapply(p.value, function(x) stars_signif(x)))) %>%
           tibble::as.tibble()
}  


## Function for ploting correlation data frames resulting from Get_Top_Relationships and Correlations_With_One()
func_dotplot_cor <- function(df){                                        # https://www.r-pkg.org/pkg/ggpubr
  dotplotcor_scale_fill <- function(...){                                # Fix colors to signif factor levels even if missing
    ggplot2:::manual_scale(                                   
      'color', 
      values = setNames(
        c("darkgreen", "green3", "lawngreen", "yellow", "red"), 
        c("***", "**", "*", ".", "ns")), 
      ...
    )
  }                                           
  
  dtoplot_theme <- 
    ggpubr::theme_pubr() +
    theme(axis.text.y = element_text(size = 10))
  
  if(!"Variable" %in% colnames(df)){                                             # in oder to work for both Get_Top_Relationships and Correlations_With_One()
  df <- 
    df %>%                                            
      unite(cor_between, c("feature_1", "feature_2"), sep = " X ")               # unite 2 columns to x name from plot
  }else df <- df %>% dplyr::rename(cor_between = Variable)                       # change Variable to x name from plot
  
  df %>%
    ggpubr::ggdotchart(x = "cor_between", y = "correlation",
                       color = "p.signif",                                       # Color by sig
                       #   palette = c("#00AFBB", "#E7B800", "#FC4E07"),         # Custom color palette
                       sorting = "descending",                                   # Sort value in descending order
                       add = "segments",                                         # Add segments from y = 0 to dots
                       add.params = list(color = "lightgray", size = 2),         # Change segment color and size
                       group = "p.signif",                                       # Order by groups
                       dot.size = 8,                                             # Large dot size
                       xlab = "",
                       rotate = TRUE,                                            # Rotate vertically
                       label = round(.$correlation, 1),                          # Add mpg values as dot labels
                       font.label = list(color = "white", size = 9, 
                                         vjust = 0.5),                           # Adjust label parameters
                       ggtheme = dtoplot_theme) +                                # ggplot2 theme
    dotplotcor_scale_fill() +                                            # Fix colors to signif factor levels even if missing
    geom_hline(yintercept = 0, linetype = 2, color = "lightgray")
}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Define Function for Pre-Post Plots, t Change and ANCOVA Post
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Func t test si boxplot simplu
func_t_box <- function(df, ind, pre_var, post_var){
  df_modif <-
    df %>%
    select(ind, pre_var, post_var) %>% 
    tidyr::drop_na() %>%
    gather(pre_var, post_var, key = "Cond", value = "value") %>% 
    mutate_at(vars(c(1, 2)), list(~as.factor(.))) %>% 
    mutate(Cond = factor(Cond, levels = c(pre_var, post_var))) 
  
  stat_comp <- ggpubr::compare_means(value ~ Cond, data = df_modif, method = "t.test", paired = TRUE)
  
  stat_comp2 <-
    df_modif %>% 
    do(tidy(t.test(.$value ~ .$Cond,
                   paired = TRUE,
                   data=.)))
  
  plot <- 
    ggpubr::ggpaired(df_modif, x = "Cond", y = "value", id = ind, 
                     color = "Cond", line.color = "gray", line.size = 0.4,
                     palette = c("#00AFBB", "#FC4E07"), legend = "none") +
      stat_summary(fun.data = mean_se,  colour = "darkred") +
      ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label.x = as.numeric(df_modif$Cond)-0.4, label.y = max(df_modif$value)+0.5) + 
      ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label = "p.signif", comparisons = list(c(pre_var, post_var)))
  
  cat(paste0("#### ", pre_var, " ", post_var, "\n", "\n"))
  print(stat_comp)
  print(stat_comp2)
  print(plot)
}


func_ancova_multibox <- function(df, ind, pre_var_c1, post_var_c1, pre_var_c2, post_var_c2){
  
  diff_score_c1 <- paste0(post_var_c1, " - ", pre_var_c1)
  diff_score_c2 <- paste0(post_var_c2, " - ", pre_var_c2)
  
  ## Plots and p-values for t tests
  df_modif <-
    df %>%
    select(ind, pre_var_c1, post_var_c1, pre_var_c2, post_var_c2) %>% 
    tidyr::drop_na() %>%
    gather(pre_var_c1, post_var_c1, pre_var_c2, post_var_c2, key = "Cond", value = "value") %>% 
    mutate_at(vars(c(1, 2)), list(~as.factor(.))) %>% 
    mutate(Cond = factor(Cond, levels = c(pre_var_c1, post_var_c1, pre_var_c2, post_var_c2))) 
  
  stat_comp <- ggpubr::compare_means(value ~ Cond, data = df_modif, method = "t.test", paired = TRUE, p.adjust.method = "holm")
  
  plot <-
    ggpubr::ggpaired(df_modif, x = "Cond", y = "value", id = ind, 
                     color = "Cond", line.color = "gray", line.size = 0.4,
                     palette = c("#00AFBB", "#FC4E07", "#00AFBB", "#FC4E07"), legend = "none") +
    stat_summary(fun.data = mean_se,  colour = "darkred") +
    ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label = "p.signif", 
                               label.y = c(max(df_modif$value) + 0.1*IQR(df_modif$value),
                                           max(df_modif$value) + 0.1*IQR(df_modif$value),
                                           seq(max(df_modif$value) + 0.3*IQR(df_modif$value), 
                                               max(df_modif$value) + 0.9*IQR(df_modif$value), length.out = 4)),  
                               comparisons = list(c(pre_var_c1, post_var_c1),
                                                  c(pre_var_c2, post_var_c2),
                                                  c(post_var_c1, pre_var_c2),
                                                  c(pre_var_c1, pre_var_c2),
                                                  c(post_var_c1, post_var_c2),
                                                  c(pre_var_c1, post_var_c2)))
  
  ## For ttestChange or ANCOVAChange - we do ttestChange (Post-Pre) here, but it isnt very important
  df_modif2 <-                                 
    df %>%
    select(ind, pre_var_c1, post_var_c1, pre_var_c2, post_var_c2) %>%
    tidyr::drop_na() 
  df_modif2[diff_score_c1] <- df_modif2[, post_var_c1] - df_modif2[, pre_var_c1]
  df_modif2[diff_score_c2] <- df_modif2[, post_var_c2] - df_modif2[, pre_var_c2]
  
  tChange <- t.test(df_modif2[, diff_score_c1], df_modif2[, diff_score_c2], paired = TRUE)
  
  ## For descriptives by 2 factors (PrePost and OglEcran)
  df_modif3 <-
    df %>%
    select(ind, pre_var_c1, post_var_c1, pre_var_c2, post_var_c2) %>%
    tidyr::drop_na() %>%
    gather(pre_var_c1, post_var_c1, pre_var_c2, post_var_c2, key = "Cond", value = "value") %>%
    mutate(PrePost = case_when(stringr::str_detect(.$"Cond", "pre|Pre") ~ "Pre",
                               stringr::str_detect(.$"Cond", "post|Post") ~ "Post",
                               TRUE ~ NA_character_),
           OglEcran = case_when(stringr::str_detect(.$"Cond", "ogl|OGL") ~ "OGL",
                              stringr::str_detect(.$"Cond", "ecran|ECRAN") ~ "ECRAN",
                              TRUE ~ NA_character_)) %>%
    mutate(PrePost = as.factor(PrePost),
           OglEcran = as.factor(OglEcran))
  
  ## For ANCOVAPost - this is what we use
  df_modif4 <-
    df_modif3 %>%
    select(-"Cond") %>%
    spread("PrePost", "value")
  
  ## Models (here we use ANCOVAPost)    # https://m-clark.github.io/docs/mixedModels/anovamixed.html#introduction
  full_ancovaPost <-                                          # this is better than using lm() and glht()
      jmv::ancova(
        formula = Post ~ Pre + OglEcran,
        data = df_modif4,
        homo = TRUE,
        ss = "3",                                                     
        postHoc = ~ OglEcran,
        postHocCorr = list("tukey"),
        effectSize = list("eta", "partEta")
      )
  
      # mod_ancovaPost <- lm(Post ~ Pre + OglEcran, data = df_modif4)            # this is a Covariate Second model
      # mod_ancovaPost_ss3 <- car::Anova(mod_ancovaPost, type = "III")         # Type III sums of squares; see Andy Fields 2012
      # postHocs <- multcomp::glht(mod_ancovaPost, linfct = multcomp::mcp(OglEcran = "Tukey"))  # differences between the adjusted means,
      # sum_postHocs <- summary(postHocs)                                               # use Tukey or Dunnett?s post hoc tests
      # conf_postHocs <- confint(postHocs)
  scatter <-                                                             # Check for homogeneity of regression slopes
    ggplot(df_modif4, aes(Pre, Post, colour = OglEcran)) +
    geom_point(aes(shape = OglEcran), size = 3) +
    geom_smooth(method = "lm", aes(fill = OglEcran), alpha = 0.1)
    
  
  ## Other Models that work for this date
  # mod_ancovaPost <- lm(post ~ pre + treat)      # exactly the same with aov(post ~ pre + treat)
  # summary(mod_ancovaPost)
  # 
  # mod_anovaRM <- aov(score ~ treat*time + Error(id), dflong)
  # summary(mod_anovaRM)
  # 
  # mod_lme <- lme4::lmer(score ~ treat*time + (1|id), data=dflong)
  # anova(lmeModel)
  
  ## Output
  print(plot)
  cat(paste0("#### ", pre_var_c1, " ", post_var_c1, " ", pre_var_c2, " ", post_var_c2, "\n", "\n"))
  
  cat("#### Descriptives")
  psych::describeBy(df_modif3[, "value"], list(df_modif3[, "PrePost"], df_modif3[, "OglEcran"]), mat = TRUE) %>% 
    as.tibble() %>%
    print()
  cat("\n")
  
  print(stat_comp)
  cat("\n")
  
  cat("#### t Change")
  tidy(tChange) %>% print()
  cat("\n")
  
  cat("#### ANCOVA Post")
  cat("\n")
  cat("##### Homogeneity test")
  print(tibble::as.tibble(full_ancovaPost$assump$homo))
  cat("##### ANCOVA output")
  print(tibble::as.tibble(full_ancovaPost$main))
  # tidy(mod_ancovaPost) %>% 
  #   mutate(p.signif = sapply(p.value, function(x) stars_signif(x))) %>% 
  #   print()
  # cat("\n")
  cat("##### Post Hoc")
  print(tibble::as.tibble(full_ancovaPost$postHoc[[1]]))
  # tidy(sum_postHocs) %>% 
  #   mutate(p.signif = sapply(p.value, function(x) stars_signif(x))) %>% 
  #   print()
  cat("\n")
  cat("##### Homogeneity of regression slopes")
  # subchunkify(plot(scatter), 5, 5)                       # THIS FUNC CRASHES PANDOC CONVERSION WITH ERROR 99 - DONT USE HERE
}

4 Analyses

4.1 Simple before-after analyses with t test

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Simple before-after analyses with t test
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## OXT has 2 outliers: 19, 43
# func_t_box(Data_merged_wide, "ID", "OXT_Pre_OGL", "OXT_Post_OGL")          # outlier in 57 sub  
# func_t_box(Data_merged_wide, "ID", "OXT_Pre_ECRAN", "OXT_Post_ECRAN")      # ns

func_t_box(Data_merged_wide_all, "ID", "OXT_Pre_OGL", "OXT_Post_OGL")    # ns

4.1.0.1 OXT_Pre_OGL OXT_Post_OGL

func_t_box(Data_merged_wide_all, "ID", "OXT_Pre_ECRAN", "OXT_Post_ECRAN")   # ns

4.1.0.2 OXT_Pre_ECRAN OXT_Post_ECRAN

func_t_box(Data_merged_wide_all, "ID", "Vas_Stres_Pre_OGL", "Vas_Stres_Post_OGL")         # sig  

4.1.0.3 Vas_Stres_Pre_OGL Vas_Stres_Post_OGL

func_t_box(Data_merged_wide_all, "ID", "Vas_Stres_Pre_ECRAN", "Vas_Stres_Post_ECRAN")     # ns

4.1.0.4 Vas_Stres_Pre_ECRAN Vas_Stres_Post_ECRAN

func_t_box(Data_merged_wide_all, "ID", "Vas_Bine_Pre_OGL", "Vas_Bine_Post_OGL")           # ns  

4.1.0.5 Vas_Bine_Pre_OGL Vas_Bine_Post_OGL

func_t_box(Data_merged_wide_all, "ID", "Vas_Bine_Pre_ECRAN", "Vas_Bine_Post_ECRAN")       # sig

4.1.0.6 Vas_Bine_Pre_ECRAN Vas_Bine_Post_ECRAN

func_t_box(Data_merged_wide_all, "ID", "IOS_Pre_OGL", "IOS_Post_OGL")                     # sig  

4.1.0.7 IOS_Pre_OGL IOS_Post_OGL

func_t_box(Data_merged_wide_all, "ID", "IOS_Pre_ECRAN", "IOS_Post_ECRAN")                 # ns

4.1.0.8 IOS_Pre_ECRAN IOS_Post_ECRAN

4.2 tChance and ANCOVAPost

## tChance and ANCOVAPost 

func_ancova_multibox(Data_merged_wide_all, "ID", "OXT_Pre_OGL", "OXT_Post_OGL", "OXT_Pre_ECRAN", "OXT_Post_ECRAN")              # ns

4.2.0.1 OXT_Pre_OGL OXT_Post_OGL OXT_Pre_ECRAN OXT_Post_ECRAN

4.2.0.2 Descriptives

NA

4.2.0.3 t Change

4.2.0.4 ANCOVA Post

4.2.0.4.1 Homogeneity test

4.2.0.4.2 ANCOVA output

4.2.0.4.3 Post Hoc

4.2.0.4.4 Homogeneity of regression slopes

func_ancova_multibox(Data_merged_wide_all,  "ID", "Vas_Stres_Pre_OGL", "Vas_Stres_Post_OGL", "Vas_Stres_Pre_ECRAN", "Vas_Stres_Post_ECRAN")    # ns

4.2.0.5 Vas_Stres_Pre_OGL Vas_Stres_Post_OGL Vas_Stres_Pre_ECRAN Vas_Stres_Post_ECRAN

4.2.0.6 Descriptives

NA

4.2.0.7 t Change

4.2.0.8 ANCOVA Post

4.2.0.8.1 Homogeneity test

4.2.0.8.2 ANCOVA output

4.2.0.8.3 Post Hoc

4.2.0.8.4 Homogeneity of regression slopes

func_ancova_multibox(Data_merged_wide_all,  "ID", "Vas_Bine_Pre_OGL", "Vas_Bine_Post_OGL", "Vas_Bine_Pre_ECRAN", "Vas_Bine_Post_ECRAN")        # ns

4.2.0.9 Vas_Bine_Pre_OGL Vas_Bine_Post_OGL Vas_Bine_Pre_ECRAN Vas_Bine_Post_ECRAN

4.2.0.10 Descriptives

NA

4.2.0.11 t Change

4.2.0.12 ANCOVA Post

4.2.0.12.1 Homogeneity test

4.2.0.12.2 ANCOVA output

4.2.0.12.3 Post Hoc

4.2.0.12.4 Homogeneity of regression slopes

func_ancova_multibox(Data_merged_wide_all,  "ID", "IOS_Pre_OGL", "IOS_Post_OGL", "IOS_Pre_ECRAN", "IOS_Post_ECRAN")                            # sig

4.2.0.13 IOS_Pre_OGL IOS_Post_OGL IOS_Pre_ECRAN IOS_Post_ECRAN

4.2.0.14 Descriptives

NA

4.2.0.15 t Change

4.2.0.16 ANCOVA Post

4.2.0.16.1 Homogeneity test

4.2.0.16.2 ANCOVA output

4.2.0.16.3 Post Hoc

4.2.0.16.4 Homogeneity of regression slopes

4.3 Simple before-after analyses with t test by Gender

# Females
cat("### Fem")

4.3.1 Fem

func_t_box(Data_fem, "ID", "OXT_Pre_OGL", "OXT_Post_OGL")    # sig

4.3.1.1 OXT_Pre_OGL OXT_Post_OGL

func_t_box(Data_fem, "ID", "OXT_Pre_ECRAN", "OXT_Post_ECRAN")   # ns

4.3.1.2 OXT_Pre_ECRAN OXT_Post_ECRAN

func_t_box(Data_fem, "ID", "Vas_Stres_Pre_OGL", "Vas_Stres_Post_OGL")         # sig 

4.3.1.3 Vas_Stres_Pre_OGL Vas_Stres_Post_OGL

func_t_box(Data_fem, "ID", "Vas_Stres_Pre_ECRAN", "Vas_Stres_Post_ECRAN")     # sig

4.3.1.4 Vas_Stres_Pre_ECRAN Vas_Stres_Post_ECRAN

func_t_box(Data_fem, "ID", "Vas_Bine_Pre_OGL", "Vas_Bine_Post_OGL")           # ns  

4.3.1.5 Vas_Bine_Pre_OGL Vas_Bine_Post_OGL

func_t_box(Data_fem, "ID", "Vas_Bine_Pre_ECRAN", "Vas_Bine_Post_ECRAN")       # ns

4.3.1.6 Vas_Bine_Pre_ECRAN Vas_Bine_Post_ECRAN

func_t_box(Data_fem, "ID", "IOS_Pre_OGL", "IOS_Post_OGL")                     # sig  

4.3.1.7 IOS_Pre_OGL IOS_Post_OGL

func_t_box(Data_fem, "ID", "IOS_Pre_ECRAN", "IOS_Post_ECRAN")                 # ns

4.3.1.8 IOS_Pre_ECRAN IOS_Post_ECRAN

# Males
cat("### Masc")

4.3.2 Masc

func_t_box(Data_masc, "ID", "OXT_Pre_OGL", "OXT_Post_OGL")    # ns

4.3.2.1 OXT_Pre_OGL OXT_Post_OGL

func_t_box(Data_masc, "ID", "OXT_Pre_ECRAN", "OXT_Post_ECRAN")   # ns

4.3.2.2 OXT_Pre_ECRAN OXT_Post_ECRAN

func_t_box(Data_masc, "ID", "Vas_Stres_Pre_OGL", "Vas_Stres_Post_OGL")         # ns  

4.3.2.3 Vas_Stres_Pre_OGL Vas_Stres_Post_OGL

func_t_box(Data_masc, "ID", "Vas_Stres_Pre_ECRAN", "Vas_Stres_Post_ECRAN")     # ns

4.3.2.4 Vas_Stres_Pre_ECRAN Vas_Stres_Post_ECRAN

func_t_box(Data_masc, "ID", "Vas_Bine_Pre_OGL", "Vas_Bine_Post_OGL")           # ns  

4.3.2.5 Vas_Bine_Pre_OGL Vas_Bine_Post_OGL

func_t_box(Data_masc, "ID", "Vas_Bine_Pre_ECRAN", "Vas_Bine_Post_ECRAN")       # ns

4.3.2.6 Vas_Bine_Pre_ECRAN Vas_Bine_Post_ECRAN

func_t_box(Data_masc, "ID", "IOS_Pre_OGL", "IOS_Post_OGL")                     # sig  

4.3.2.7 IOS_Pre_OGL IOS_Post_OGL

func_t_box(Data_masc, "ID", "IOS_Pre_ECRAN", "IOS_Post_ECRAN")                 # ns

4.3.2.8 IOS_Pre_ECRAN IOS_Post_ECRAN

4.4 tChance and ANCOVAPost by Gender

# Females
cat("### Fem")

4.4.1 Fem

Data_fem %>%
  func_ancova_multibox("ID", "OXT_Pre_OGL", "OXT_Post_OGL", "OXT_Pre_ECRAN", "OXT_Post_ECRAN")                                    # ns (PROBLEM)

4.4.1.1 OXT_Pre_OGL OXT_Post_OGL OXT_Pre_ECRAN OXT_Post_ECRAN

4.4.1.2 Descriptives

NA

4.4.1.3 t Change

4.4.1.4 ANCOVA Post

4.4.1.4.1 Homogeneity test

4.4.1.4.2 ANCOVA output

4.4.1.4.3 Post Hoc

4.4.1.4.4 Homogeneity of regression slopes

func_ancova_multibox(Data_fem, "ID", "Vas_Stres_Pre_OGL", "Vas_Stres_Post_OGL", "Vas_Stres_Pre_ECRAN", "Vas_Stres_Post_ECRAN")    # ns

4.4.1.5 Vas_Stres_Pre_OGL Vas_Stres_Post_OGL Vas_Stres_Pre_ECRAN Vas_Stres_Post_ECRAN

4.4.1.6 Descriptives

NA

4.4.1.7 t Change

4.4.1.8 ANCOVA Post

4.4.1.8.1 Homogeneity test

4.4.1.8.2 ANCOVA output

4.4.1.8.3 Post Hoc

4.4.1.8.4 Homogeneity of regression slopes

func_ancova_multibox(Data_fem, "ID", "Vas_Bine_Pre_OGL", "Vas_Bine_Post_OGL", "Vas_Bine_Pre_ECRAN", "Vas_Bine_Post_ECRAN")        # ns

4.4.1.9 Vas_Bine_Pre_OGL Vas_Bine_Post_OGL Vas_Bine_Pre_ECRAN Vas_Bine_Post_ECRAN

4.4.1.10 Descriptives

NA

4.4.1.11 t Change

4.4.1.12 ANCOVA Post

4.4.1.12.1 Homogeneity test

4.4.1.12.2 ANCOVA output

4.4.1.12.3 Post Hoc

4.4.1.12.4 Homogeneity of regression slopes

func_ancova_multibox(Data_fem, "ID", "IOS_Pre_OGL", "IOS_Post_OGL", "IOS_Pre_ECRAN", "IOS_Post_ECRAN")                            # sig

4.4.1.13 IOS_Pre_OGL IOS_Post_OGL IOS_Pre_ECRAN IOS_Post_ECRAN

4.4.1.14 Descriptives

NA

4.4.1.15 t Change

4.4.1.16 ANCOVA Post

4.4.1.16.1 Homogeneity test

4.4.1.16.2 ANCOVA output

4.4.1.16.3 Post Hoc

4.4.1.16.4 Homogeneity of regression slopes

# Males
cat("### Masc")

4.4.2 Masc

Data_masc %>%
  func_ancova_multibox("ID", "OXT_Pre_OGL", "OXT_Post_OGL", "OXT_Pre_ECRAN", "OXT_Post_ECRAN")                                     # ns

4.4.2.1 OXT_Pre_OGL OXT_Post_OGL OXT_Pre_ECRAN OXT_Post_ECRAN

4.4.2.2 Descriptives

NA

4.4.2.3 t Change

4.4.2.4 ANCOVA Post

4.4.2.4.1 Homogeneity test

4.4.2.4.2 ANCOVA output

4.4.2.4.3 Post Hoc

4.4.2.4.4 Homogeneity of regression slopes

func_ancova_multibox(Data_masc, "ID", "Vas_Stres_Pre_OGL", "Vas_Stres_Post_OGL", "Vas_Stres_Pre_ECRAN", "Vas_Stres_Post_ECRAN")    # ns

4.4.2.5 Vas_Stres_Pre_OGL Vas_Stres_Post_OGL Vas_Stres_Pre_ECRAN Vas_Stres_Post_ECRAN

4.4.2.6 Descriptives

NA

4.4.2.7 t Change

4.4.2.8 ANCOVA Post

4.4.2.8.1 Homogeneity test

4.4.2.8.2 ANCOVA output

4.4.2.8.3 Post Hoc

4.4.2.8.4 Homogeneity of regression slopes

func_ancova_multibox(Data_masc, "ID", "Vas_Bine_Pre_OGL", "Vas_Bine_Post_OGL", "Vas_Bine_Pre_ECRAN", "Vas_Bine_Post_ECRAN")        # ns

4.4.2.9 Vas_Bine_Pre_OGL Vas_Bine_Post_OGL Vas_Bine_Pre_ECRAN Vas_Bine_Post_ECRAN

4.4.2.10 Descriptives

NA

4.4.2.11 t Change

4.4.2.12 ANCOVA Post

4.4.2.12.1 Homogeneity test

4.4.2.12.2 ANCOVA output

4.4.2.12.3 Post Hoc

4.4.2.12.4 Homogeneity of regression slopes

func_ancova_multibox(Data_masc, "ID", "IOS_Pre_OGL", "IOS_Post_OGL", "IOS_Pre_ECRAN", "IOS_Post_ECRAN")                            # sig

4.4.2.13 IOS_Pre_OGL IOS_Post_OGL IOS_Pre_ECRAN IOS_Post_ECRAN

4.4.2.14 Descriptives

NA

4.4.2.15 t Change

4.4.2.16 ANCOVA Post

4.4.2.16.1 Homogeneity test

4.4.2.16.2 ANCOVA output

4.4.2.16.3 Post Hoc

4.4.2.16.4 Homogeneity of regression slopes

4.5 Not considering Cond (OglEcran) - Overall Pre-Post

# Data
vars_OglEcran <- grep("OGL|ECRAN", colnames(Data_merged_wide_all))

Data_notOglEcran <-
  Data_merged_wide_all %>%
  gather(variable, value, vars_OglEcran, -c(ID)) %>%                                      # gather Data and Nr_zi as these are OglEcran level
  tidyr::separate(variable,  c("variable", "OglEcran"), "_(?=[^_]+$)") %>%                # split only on last "_"
  spread(variable, value) %>% 
  rename_at(vars(vars_OglEcran), list(~stringr::str_remove_all(., c("_OGL|_ECRAN")))) %>%   
  mutate_at(vars(c(45:46, 48:53)), list(~as.numeric(as.character(.)))) %>%
  mutate(Diff_OXT = OXT_Post - OXT_Pre,
         Diff_Vas_Stres = Vas_Stres_Post - Vas_Stres_Pre,
         Diff_Vas_Bine = Vas_Bine_Post - Vas_Bine_Pre
  ) %>%
  arrange(ID)
  

# Simple before-after analyses with t test Not considering Cond (OglEcran)
Data_notOglEcran %>%
  filter(!ID %in% c(19, 43)) %>%
  func_t_box("ID", "OXT_Pre", "OXT_Post")    # JUST SHY OF SIG

4.5.0.1 OXT_Pre OXT_Post

func_t_box(Data_notOglEcran, "ID", "Vas_Stres_Pre", "Vas_Stres_Post")         # sig  

4.5.0.2 Vas_Stres_Pre Vas_Stres_Post

func_t_box(Data_notOglEcran, "ID", "Vas_Bine_Pre", "Vas_Bine_Post")           # sig  

4.5.0.3 Vas_Bine_Pre Vas_Bine_Post

func_t_box(Data_notOglEcran, "ID", "IOS_Pre", "IOS_Post")                     # sig  

4.5.0.4 IOS_Pre IOS_Post

# GGpairs
fn_ggpair_reg <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    geom_point() + 
    geom_smooth(method=loess, fill="red", color="red", ...) +
    geom_smooth(method=lm, fill="blue", color="blue", ...)
  p
}

GGally::ggpairs(data = Data_notOglEcran %>% filter(!ID %in% c(19, 43)),
                columns = 54:56, 
                lower = list(continuous = fn_ggpair_reg))

4.6 Difference scores conditioned on Anxiety ?

plot_columns <- match( c("StaiSbrut", "StaiTbrut",
                         "Diff_OXT_OGL", "Diff_OXT_ECRAN",
                         "Diff_Vas_Stres_OGL", "Diff_Vas_Stres_ECRAN",
                         "Diff_Vas_Bine_OGL", "Diff_Vas_Bine_ECRAN"), 
                          names(Data_difscores))

# GGpairs
fn_ggpair_reg <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    geom_point() + 
    geom_smooth(method=loess, fill="red", color="red", ...) +
    geom_smooth(method=lm, fill="blue", color="blue", ...)
  p
}

GGally::ggpairs(data = Data_difscores,
                columns = plot_columns, 
                lower = list(continuous = fn_ggpair_reg))

# GGpairs By Gen
GGally::ggpairs(data = Data_difscores,
                columns = plot_columns, 
                title = "correlation matrix",               
                mapping = aes(colour = Gen), 
                lower = list(
                continuous = "smooth",
                combo = "facetdensity",
                mapping = aes(color = Gen)),
                upper = list(continuous = wrap("cor", size = 3, hjust=0.8)))

# Coplot Anx 
Data_difscores %>%
  filter(!ID %in% c(19, 43)) %>%
    coplot(Diff_OXT_OGL ~ Diff_Vas_Stres_OGL | StaiTbrut,
           data = .,
           rows = 1,
           panel = function(x, y, ...) {
           panel.smooth(x, y, span = .8, iter = 5,...)
           abline(lm(y ~ x), col = "blue")})

Missing rows: 17, 26, 34, 43

Data_difscores %>%
  filter(!ID %in% c(19, 43)) %>%
    coplot(Diff_OXT_ECRAN ~ Diff_Vas_Stres_ECRAN | StaiTbrut,
           data = .,
           rows = 1,
           panel = function(x, y, ...) {
           panel.smooth(x, y, span = .8, iter = 5,...)
           abline(lm(y ~ x), col = "blue")})

Missing rows: 17, 26, 34, 43

Data_difscores %>%
  filter(!ID %in% c(19, 43)) %>%
  coplot(Diff_OXT_OGL ~ Diff_Vas_Stres_OGL | StaiSbrut,
         data = .,
         rows = 1,
         panel = function(x, y, ...) {
           panel.smooth(x, y, span = .8, iter = 5,...)
           abline(lm(y ~ x), col = "blue")})

Missing rows: 17, 26, 34, 43

Data_difscores %>%
  filter(!ID %in% c(19, 43)) %>%
  coplot(Diff_OXT_ECRAN ~ Diff_Vas_Stres_ECRAN | StaiSbrut,
         data = .,
         rows = 1,
         panel = function(x, y, ...) {
           panel.smooth(x, y, span = .8, iter = 5,...)
           abline(lm(y ~ x), col = "blue")})

Missing rows: 17, 26, 34, 43

4.7 Difference scores conditioned on Extraversion ?

plot_columns <- match(c("E",
                        "Diff_OXT_OGL", "Diff_OXT_ECRAN",
                        "Diff_Vas_Stres_OGL", "Diff_Vas_Stres_ECRAN",
                        "Diff_Vas_Bine_OGL", "Diff_Vas_Bine_ECRAN"), 
                        names(Data_difscores))

# GGpairs
fn_ggpair_reg <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    geom_point() + 
    geom_smooth(method=loess, fill="red", color="red", ...) +
    geom_smooth(method=lm, fill="blue", color="blue", ...)
  p
}

GGally::ggpairs(data = Data_difscores,
                columns = plot_columns, 
                lower = list(continuous = fn_ggpair_reg))

## Define Function
func_moderation <- function(Data, dep, mod, pred){
  moderation <- 
    Data %>%
    medmod::mod(., dep = dep, mod = mod, pred = pred,
                ci = TRUE, estMethod = 'standard', test = TRUE, simpleSlopeEst = FALSE, simpleSlopePlot = TRUE)
  
  cat(paste("<b> Moderation: ", "Dep = ", dep, "Pred = ", pred, "Mod = ", mod, "</b>"))
  moderation$mod %>% 
    knitr::kable(caption = "Moderation", digits = 3) %>%
    print()
  moderation$simpleSlope$plot %>%
    print()
}


Data_merged_wide_all %>%
  filter(!ID %in% c(19, 43)) %>%
  func_moderation(Data = ., dep = "OXT_Post_OGL", mod = "E", pred = "OXT_Pre_OGL")

Data_merged_wide_all %>%
  filter(!ID %in% c(19, 43)) %>%
  func_moderation(Data = ., dep = "OXT_Post_ECRAN", mod = "E", pred = "OXT_Pre_ECRAN")


cat("### ANCOVA - Extraversion")
df_ancova_e <-
  Data_merged_wide_all %>%
  select("ID", "OXT_Pre_OGL", "OXT_Post_OGL", "OXT_Pre_ECRAN", "OXT_Post_ECRAN", "E") %>%
  tidyr::drop_na() %>%
  gather("OXT_Pre_OGL", "OXT_Post_OGL", "OXT_Pre_ECRAN", "OXT_Post_ECRAN", key = "Cond", value = "value") %>%
  mutate(PrePost = case_when(stringr::str_detect(.$"Cond", "Pre") ~ "Pre",
                             stringr::str_detect(.$"Cond", "Post") ~ "Post",
                             TRUE ~ NA_character_),
         OglEcran = case_when(stringr::str_detect(.$"Cond", "OGL") ~ "OGL",
                              stringr::str_detect(.$"Cond", "ECRAN") ~ "ECRAN",
                              TRUE ~ NA_character_)) %>%
  mutate(PrePost = as.factor(PrePost),
         OglEcran = as.factor(OglEcran)) %>%
  select(-"Cond") %>%
  spread("PrePost", "value")


jmv::ancova(
    formula = Post ~ Pre + OglEcran + E,
    data = df_ancova_e,
    homo = TRUE,
    ss = "3",
    postHoc = ~ OglEcran,
    postHocCorr = list("tukey"),
    effectSize = list("eta", "partEta")
  )

## Define Function
func_moderation <- function(Data, dep, mod, pred){
  moderation <- 
    Data %>%
      medmod::mod(., dep = dep, mod = mod, pred = pred,
                  ci = TRUE, estMethod = 'standard', test = TRUE, simpleSlopeEst = FALSE, simpleSlopePlot = TRUE)
  
  cat(paste("<b> Moderation: ", "Dep = ", dep, "Pred = ", pred, "Mod = ", mod, "</b>"))
  moderation$mod %>% 
    knitr::kable(caption = "Moderation", digits = 3) %>%
    print()
  moderation$simpleSlope$plot %>%
    print()
}

## Apply Function
Data_merged_wide_all %>%
  filter(!ID %in% c(19, 43)) %>%
  func_moderation(Data = ., dep = "OXT_Post_OGL", mod = "StaiTbrut", pred = "OXT_Pre_OGL")

Data_merged_wide_all %>%
  filter(!ID %in% c(19, 43)) %>%
  func_moderation(Data = ., dep = "OXT_Post_ECRAN", mod = "StaiTbrut", pred = "OXT_Pre_ECRAN")





# bla <-
#   Data_merged %>%
#   mutate(Conditie = as.numeric(as.factor(Conditie)))
# 
# psych::mediate(data = bla, Oxitocina_post ~ Oxitocina_pre + Conditie + Vas_rel_global)   # moderation with covariate and diagram

4.8 Mixed Design ANOVA

# Data
Data_mixedanova <-
  Data_merged_long_all %>%
  mutate(PrePost = case_when(stringr::str_detect(.$"PrePost", "Pre") ~ as.integer(0),
                             stringr::str_detect(.$"PrePost", "Post") ~ as.integer(1),
                             TRUE ~ NA_integer_)) %>%
  mutate(OglEcran = case_when(stringr::str_detect(.$"OglEcran", "ECRAN") ~ as.integer(0),
                           stringr::str_detect(.$"OglEcran", "OGL") ~ as.integer(1),
                           TRUE ~ NA_integer_)) %>%
  mutate(StaiTbrut_dih = sjmisc::dicho(.$"StaiTbrut", 
                                       dich.by = "median", as.num = FALSE,                                # similar to Hmisc::cut2, but uses median
                                       val.labels = c("lower", "higher"), append = FALSE)) %>%
  select(ID, OglEcran, PrePost, OXT, StaiTbrut_dih) %>%
  drop_na()
           
    
# Mixed Design Anova  # http://www.cookbook-r.com/Statistical_analysis/ANOVA/
aov_mixed1 <- aov(OXT ~ OglEcran * PrePost * StaiTbrut_dih + Error(ID/PrePost), data = Data_mixedanova) 
broom::tidy(aov_mixed1)
# summary(aov_mixed1)   
# model.tables(aov_mixed1, "means")   # is data balanced?




# library(jtools)
# mod1 <- lm(OXT_Post_OGL ~ OXT_Pre_OGL * StaiTbrut, data = Data_mixedanova)
# jtools::summ(mod1)
# jtools::interact_plot(mod1, pred = OXT_Pre_OG, modx = StaiTbrut)
# 
# library(interactions)

–>

4.9 Correlations between Diff_OXT_OGL with other variables

## Correlations between Diffrence Scores with other variables
Data_difscores[,-c(1:2, 5:24)] %>%
  Correlations_With_One(., variable = "Diff_OXT_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05)   # nothing

4.9.0.1 Correlations with Diff_OXT_OGL

####
Data_difscores[,-c(1:2, 5:24)] %>%
   Correlations_With_One(., variable = "Diff_Vas_Stres_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05) # C,C5

4.9.0.2 Correlations with Diff_Vas_Stres_OGL

###
Data_difscores[,-c(1:2, 5:24)] %>%
   Correlations_With_One(., variable = "Diff_Vas_Bine_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05)  # BDI

4.9.0.3 Correlations with Diff_Vas_Bine_OGL

Data_difscores[,-c(1:2, 5:24)] %>%
   Correlations_With_One(., variable = "Diff_Vas_Bine_ECRAN", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05) # O2,O3, StaiT

4.9.0.4 Correlations with Diff_Vas_Bine_ECRAN

lm(Vas_Stres_Post_OGL ~ Vas_Stres_Pre_OGL + StaiTbrut, data = Data_merged_wide_all) %>% summary()      # sig

Call: lm(formula = Vas_Stres_Post_OGL ~ Vas_Stres_Pre_OGL + StaiTbrut, data = Data_merged_wide_all)

Residuals: Min 1Q Median 3Q Max -18.513 -4.540 1.046 3.503 20.137

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.31395 5.38906 -1.728 0.0905 .
Vas_Stres_Pre_OGL 0.73616 0.05491 13.406 <0.0000000000000002 ** StaiTbrut 0.26002 0.12212 2.129 0.0385
— Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.039 on 47 degrees of freedom Multiple R-squared: 0.8065, Adjusted R-squared: 0.7983 F-statistic: 97.97 on 2 and 47 DF, p-value: < 0.00000000000000022

lm(Vas_Stres_Post_ECRAN ~ Vas_Stres_Pre_ECRAN + StaiTbrut, data = Data_merged_wide_all) %>% summary()  # ns

Call: lm(formula = Vas_Stres_Post_ECRAN ~ Vas_Stres_Pre_ECRAN + StaiTbrut, data = Data_merged_wide_all)

Residuals: Min 1Q Median 3Q Max -23.0762 -4.4267 -0.5844 4.6774 19.2660

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.34215 6.03290 1.217 0.230
Vas_Stres_Pre_ECRAN 0.83186 0.05658 14.702 <0.0000000000000002 *** StaiTbrut -0.13710 0.14472 -0.947 0.348
— Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.858 on 46 degrees of freedom (1 observation deleted due to missingness) Multiple R-squared: 0.8373, Adjusted R-squared: 0.8302 F-statistic: 118.4 on 2 and 46 DF, p-value: < 0.00000000000000022

lm(Vas_Bine_Post_OGL ~ Vas_Bine_Pre_OGL + ScorBDI, data = Data_merged_wide_all) %>% summary()        # sig

Call: lm(formula = Vas_Bine_Post_OGL ~ Vas_Bine_Pre_OGL + ScorBDI, data = Data_merged_wide_all)

Residuals: Min 1Q Median 3Q Max -31.928 -7.000 0.981 6.480 31.480

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.69505 6.62844 2.368 0.0221 *
Vas_Bine_Pre_OGL 0.89920 0.08809 10.208 0.000000000000164 ** ScorBDI -0.62074 0.26139 -2.375 0.0217
— Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.54 on 47 degrees of freedom Multiple R-squared: 0.7015, Adjusted R-squared: 0.6888 F-statistic: 55.24 on 2 and 47 DF, p-value: 0.0000000000004571

lm(Vas_Bine_Post_ECRAN ~ Vas_Bine_Pre_ECRAN + ScorBDI, data = Data_merged_wide_all) %>% summary()    # ns

Call: lm(formula = Vas_Bine_Post_ECRAN ~ Vas_Bine_Pre_ECRAN + ScorBDI, data = Data_merged_wide_all)

Residuals: Min 1Q Median 3Q Max -21.9414 -4.4338 -0.7211 5.6113 23.7907

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.67270 5.16231 2.067 0.0443 *
Vas_Bine_Pre_ECRAN 0.90116 0.06273 14.366 <0.0000000000000002 *** ScorBDI -0.10869 0.20274 -0.536 0.5945
— Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.741 on 46 degrees of freedom (1 observation deleted due to missingness) Multiple R-squared: 0.8259, Adjusted R-squared: 0.8183 F-statistic: 109.1 on 2 and 46 DF, p-value: < 0.00000000000000022

4.10 Correlations between Diff_OXT_OGL with other variables - by Gender

## Correlations between Diffrence Scores with other variables - Females
cat("### Females")

4.10.1 Females

Data_difscores[,-c(1:2, 5:24)] %>%
  filter(Gen == "f") %>%
  Correlations_With_One(., variable = "Diff_OXT_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05)   # nothing

4.10.1.1 Correlations with Diff_OXT_OGL

####
Data_difscores[,-c(1:2, 5:24)]%>%
  filter(Gen == "f") %>%
   Correlations_With_One(., variable = "Diff_Vas_Stres_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05) # A1

4.10.1.2 Correlations with Diff_Vas_Stres_OGL

###
Data_difscores[,-c(1:2, 5:24)] %>%
  filter(Gen == "f") %>%
  Correlations_With_One(., variable = "Diff_Vas_Bine_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05)  # BDI, N2

4.10.1.3 Correlations with Diff_Vas_Bine_OGL

Data_difscores[,-c(1:2, 5:24)] %>%
  filter(Gen == "f") %>%
  Correlations_With_One(., variable = "Diff_Vas_Bine_ECRAN", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05) # StaiT

4.10.1.4 Correlations with Diff_Vas_Bine_ECRAN

## Correlations between Diffrence Scores with other variables - Males
cat("### Males")

4.10.2 Males

Data_difscores[,-c(1:2, 5:24)] %>%
  filter(Gen == "m") %>%
  Correlations_With_One(., variable = "Diff_OXT_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05)   # O6,A6

4.10.2.1 Correlations with Diff_OXT_OGL

####
Data_difscores[,-c(1:2, 5:24)] %>%
  filter(Gen == "m") %>%
  Correlations_With_One(., variable = "Diff_Vas_Stres_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05) # C2-Ordine!!!

4.10.2.2 Correlations with Diff_Vas_Stres_OGL

###
Data_difscores[,-c(1:2, 5:24)] %>%
  filter(Gen == "m") %>%
  Correlations_With_One(., variable = "Diff_Vas_Bine_OGL", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05)  # nothing

4.10.2.3 Correlations with Diff_Vas_Bine_OGL

Data_difscores[,-c(1:2, 5:24)] %>%
  filter(Gen == "m") %>%
  Correlations_With_One(., variable = "Diff_Vas_Bine_ECRAN", correlation_abs_threshold = 0.2, pvalue_threshold = 0.05) # C2-Ordine,O2

4.10.2.4 Correlations with Diff_Vas_Bine_ECRAN

5 within-participant mediation

library(JSmediation)

df_cordif_fem <- 
  Data_difscores %>%
  filter(Gen == "f")

cor.test(df_cordif_fem$Diff_OXT_OGL, df_cordif_fem$Diff_Vas_Stres_OGL)  # ns

Pearson's product-moment correlation

data: df_cordif_fem\(Diff_OXT_OGL and df_cordif_fem\)Diff_Vas_Stres_OGL t = 0.21086, df = 22, p-value = 0.8349 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.3651008 0.4403297 sample estimates: cor 0.04491029

summary(lm(Diff_OXT_OGL ~ Diff_Vas_Stres_OGL, data = df_cordif_fem))   # ns

Call: lm(formula = Diff_OXT_OGL ~ Diff_Vas_Stres_OGL, data = df_cordif_fem)

Residuals: Min 1Q Median 3Q Max -0.62386 -0.15832 -0.07481 0.09697 1.18405

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.176485 0.089349 1.975 0.0609 . Diff_Vas_Stres_OGL 0.001570 0.007445 0.211 0.8349
— Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3689 on 22 degrees of freedom (3 observations deleted due to missingness) Multiple R-squared: 0.002017, Adjusted R-squared: -0.04335 F-statistic: 0.04446 on 1 and 22 DF, p-value: 0.8349

withinmed_OGL <- 
  mdt_within_wide(data = Data_fem, 
                  DV_A = Vas_Stres_Pre_OGL, DV_B = Vas_Stres_Post_OGL, 
                  M_A = OXT_Pre_OGL, M_B = OXT_Post_OGL)
add_index(withinmed_OGL)

6 Test of mediation (within-participant mediation)

Variables:

• DV difference: Vas_Stres_Pre_OGL - Vas_Stres_Post_OGL
• M difference: OXT_Pre_OGL - OXT_Post_OGL

Paths:

==== ============== ===== ====================== Path Point estimate SE APA
==== ============== ===== ====================== a -0.166 0.074 t(23) = 2.26, p = .034 b 3.656 6.764 t(21) = 0.54, p = .595 c 6.185 2.137 t(26) = 2.89, p = .008 c’ 7.067 2.444 t(21) = 2.89, p = .009 ==== ============== ===== ======================

Indirect effect index:

• type: Within-participant indirect effect
• point estimate: -0.608
• confidence interval:
  • method: Monte Carlo (5000 iterations)
  • level: 0.05
  • CI: [-3.37; 1.77]

Fitted models:

• 1 -> DV_diff
• 1 -> M_diff
• 1 + M_diff + M_mean -> DV_diff

withinmed_ECRAN <- 
  mdt_within_wide(data = Data_fem, 
                  DV_A = Vas_Stres_Pre_ECRAN, DV_B = Vas_Stres_Post_ECRAN, 
                  M_A = OXT_Pre_ECRAN, M_B = OXT_Post_ECRAN)
add_index(withinmed_ECRAN)

7 Test of mediation (within-participant mediation)

Variables:

• DV difference: Vas_Stres_Pre_ECRAN - Vas_Stres_Post_ECRAN
• M difference: OXT_Pre_ECRAN - OXT_Post_ECRAN

Paths:

==== ============== ===== ====================== Path Point estimate SE APA
==== ============== ===== ====================== a -0.058 0.062 t(23) = 0.93, p = .360 b 4.307 5.702 t(21) = 0.76, p = .458 c 4.923 1.988 t(25) = 2.48, p = .020 c’ 3.417 1.731 t(21) = 1.97, p = .062 ==== ============== ===== ======================

Indirect effect index:

• type: Within-participant indirect effect
• point estimate: -0.25
• confidence interval:
  • method: Monte Carlo (5000 iterations)
  • level: 0.05
  • CI: [-1.68; 0.631]

Fitted models:

• 1 -> DV_diff
• 1 -> M_diff
• 1 + M_diff + M_mean -> DV_diff

# mdt_within_wide(data = Data_fem,                                             # ns
#                 DV_A = Vas_Bine_Pre_OGL, DV_B = Vas_Bine_Post_OGL,
#                 M_A = OXT_Pre_OGL, M_B = OXT_Post_OGL)

# mdt_within_wide(data = Data_fem,                                             # ns
#                 DV_A = IOS_Pre_OGL, DV_B = IOS_Post_OGL,
#                 M_A = OXT_Pre_OGL, M_B = OXT_Post_OGL)

7.1 Correlations with OxT

## Mining Correlations with Oxy

# df_OxyAll_cor <-
#   Get_Top_Relationships(Data[,-c(1:7)], correlation_abs_threshold = 0.2, pvalue_threshold = 0.05) %>%
#     dplyr::as_tibble() %>%
#     dplyr::filter_all(any_vars(grepl("Ox", .)))                           # only Oxy, but both Pre and Post Oxy
# 
# df_OxyAll_cor %>%              
#     print(n = Inf)                                     
# df_OxyAll_cor %>%
#   func_dotplot_cor()
#  
# # Correlations only on Pre measures with Oxy -- nothing
# df_OxyPre_cor <- 
#   Get_Top_Relationships(Data[,-c(1:7)], correlation_abs_threshold = 0.2, pvalue_threshold = 0.1) %>%
#     dplyr::as_tibble() %>%
#     filter_at(vars(feature_1, feature_2), all_vars(grepl("pre|Pre", .))) %>%
#     dplyr::filter_all(any_vars(grepl("Ox", .))) 

8 Session Info

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8.1 x64 (build 9600)

Matrix products: default

locale:
[1] LC_COLLATE=Romanian_Romania.1250  LC_CTYPE=Romanian_Romania.1250    LC_MONETARY=Romanian_Romania.1250 LC_NUMERIC=C                     
[5] LC_TIME=Romanian_Romania.1250    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] JSmediation_0.1.0  GGally_1.4.0       rio_0.5.16         plyr_1.8.4         summarytools_0.8.8 DT_0.5             ggpubr_0.2        
 [8] magrittr_1.5       broom_0.5.2        papaja_0.1.0.9842  psych_1.8.12       forcats_0.4.0      stringr_1.4.0      dplyr_0.8.3       
[15] purrr_0.3.2        readr_1.3.1        tidyr_1.0.0        tibble_2.1.3       ggplot2_3.2.1      tidyverse_1.2.1    pacman_0.5.1      

loaded via a namespace (and not attached):
 [1] TH.data_1.0-9      jmvcore_1.0.8      colorspace_1.4-1   ggsignif_0.4.0     rjson_0.2.20       pryr_0.1.4         ellipsis_0.3.0    
 [8] estimability_1.3   base64enc_0.1-3    rstudioapi_0.8     mvtnorm_1.0-11     lubridate_1.7.4    xml2_1.2.0         codetools_0.2-16  
[15] splines_3.6.1      mnormt_1.5-5       knitr_1.25         zeallot_0.1.0      jsonlite_1.6       shiny_1.2.0        compiler_3.6.1    
[22] httr_1.4.0         emmeans_1.4.1      backports_1.1.4    assertthat_0.2.1   Matrix_1.2-17      lazyeval_0.2.2     cli_1.1.0         
[29] later_0.7.5        prettyunits_1.0.2  htmltools_0.3.6    tools_3.6.1        coda_0.19-2        gtable_0.3.0       glue_1.3.1        
[36] reshape2_1.4.3     Rcpp_1.0.2         carData_3.0-2      cellranger_1.1.0   vctrs_0.2.0        nlme_3.1-140       crosstalk_1.0.0   
[43] xfun_0.9           openxlsx_4.1.0     rvest_0.3.2        mime_0.7           lifecycle_0.1.0    MASS_7.3-51.4      zoo_1.8-4         
[50] scales_1.0.0       hms_0.5.1          promises_1.0.1     parallel_3.6.1     sandwich_2.5-0     RColorBrewer_1.1-2 yaml_2.2.0        
[57] curl_3.2           pander_0.6.3       reshape_0.8.8      stringi_1.4.3      highr_0.8          zip_1.0.0          rlang_0.4.0       
[64] pkgconfig_2.0.3    matrixStats_0.54.0 bitops_1.0-6       evaluate_0.14      lattice_0.20-38    rapportools_1.0    htmlwidgets_1.3   
[71] labeling_0.3       tidystats_0.3      tidyselect_0.2.5   R6_2.4.0           generics_0.0.2     multcomp_1.4-8     pillar_1.4.2      
[78] haven_2.1.1        foreign_0.8-71     withr_2.1.2        survival_2.44-1.1  abind_1.4-5        RCurl_1.95-4.11    modelr_0.1.5      
[85] crayon_1.3.4       car_3.0-2          rmarkdown_1.17     progress_1.2.0     grid_3.6.1         readxl_1.1.0       data.table_1.11.8 
[92] digest_0.6.21      xtable_1.8-4       httpuv_1.4.5       munsell_0.5.0      jmv_1.0.8         

 

---

A work by Claudiu Papasteri

claudiu.papasteri@gmail.com